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ABSTRACT 

We present a measurement of the angular power spectrum of the cosmic far-infrared background 
(CFIRB) anisotropies in one of the extragalactic fields of the Herschel Astrophysical Terahertz Large 
Area Survey (H-ATLAS) at 250, 350 and 500 /im bands. Consistent with recent measurements of 
the CFIRB power spectrum in Herschel-SPIRE maps, we confirm the existence of a clear one-halo 
term of galaxy clustering on arcminute angular scales with large-scale two-halo term of clustering 
at 30 arcminutes to angular scales of a few degrees. The power spectrum at the largest angular 
scales, especially at 250 ^m, is contaminated by the Galactic cirrus. The angular power spectrum 
is modeled using a conditional luminosity function approach to describe the spatial distribution of 
unresolved galaxies that make up the bulk of the CFIRB. Integrating over the dusty galaxy population 
responsible for the background anisotropics, we find that the cosmic abundance of dust, relative to the 
critical density, to be between ildust = 10~^ and 8 x 10~^ in the redshift range z ~ — 3. This dust 
abundance is consistent with estimates of the dust content in the Universe using quasar reddening 
and magnification measurements in the SDSS. 

Subject headings: cosmology: observations — submillimeter: galaxies — infrared: galaxies — galaxies: 
evolution — cosmology: large-scale structure of Universe 



I. INTRODUCTION 

While the total intensity of the cosmic far infrared 
background (CFIR B) is known from absolute photom- 
etry measuremen ts (jPuget et al.lll996l : iFixen et al.lll998l : 
iDwek et al.l[l99l . we still lack a complete knowledge of 
the sources, in the form of dusty star- forming galaxies, 
that make up the background. Limited by aperture sizes 
and the resulting source confusion noise (jNguven et al.l 
|2010( ). existing deep surveys with the Herschel Space Ob- 
servatory Q (Pilbratt et al. 2010) and ground-based sub- 
mm and mm- wave instruments resolve anywhere between 
5 and 15% of the background into ind ividual galaxies 
(ICoppin et al. 2006; Scott ct al. 2010; Ol iver et al.ll2010l: 
iClementset al.ll2010t iBerta et al.ii2011l ). Anisotropies of 
the CFIRB, or the spatial fluctuations of the background 
intensity, provide additional statistical information on 

^ Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with 
important participation from NASA. 



the fainter galaxies, especially those that make up the 
bulk of the background. 

While the fainter galaxies are individually undetected, 
due to gravitational growth and evolution in the large- 
scale structure t hese galaxies are exp ected to be clustered 
filMaddox et al.ll2010l : iHickox et aLl l2010': 'Kampen et al' 



2013). In the ansatz of the halo model (Cooray & Shot 
20021 ) such clustering of galaxies captures certain prop- 
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erties of the dark matter halos in which galaxies are 
found and the statistics of how those galaxies occupy 
the dark matter halos. The resulting anisotropics of 
the CFIRB are then a reflection of the spatial clus- 
tering of galaxies and their infrared luminosity. These 
CFIRB anisotropies, are best studied from the angular 
power spectrum of the background infrared light. Sep- 
aratel y, statistics such a s the probability of deflection, 
P{D) (jGl enn et al.ll2010f ). probe the variance and higher 
order cumulant statistics of the intensity variations at 
the beam scale. 
While early attempts to measure the angular power 



2 



spectrum of th e CFIRB resulted in low signal-to-nois e 
measurements (|Lagache et alj [2007t IViero et all [20091 ) . 
a first clear detection of the CFIRB power spectrum 
with Herschel-SFIRE (Griffin et al. 2010) maps be- 
tween 30 arcse conds and 30 arcmin ute angular scales 
was reported in lAmblard et alJ ()2011l ) . Those first mea- 
surements also confirmed the interpretation that galax- 
ies at the peak epoch of star formation in the Universe 
at redshifts of 1 to 3 trace the underlying dark matter 
halo distribution. Since then, additional measurements 
of the CFIRB power spect rum have come from Planck 
(jPlanck collaborationll2011| ) an d with additional SPIRE 
maps from the HerMES survey (jViero et al.ll2012f ). With 
multiple fields spanning up to 20 deg^, recent HerMES 
CFIRB power spectra probe angular scales of about 
30" to 2°. The halo model interpretation of the Her- 
MES spectra suggest that the halo mass scale for peak 
star-formation activity is logMpoak/A^o ^ 13.9 ± 0.6 
and the minimum halo mass to host dusty galaxies is 
log Afmi„/M0 - 10.8 ± 0.6. 

The angular power spectrum of CFIRB, in principle, 
captures the spatial distribution of the background inten- 
sity, regardless of whether the emission is from individual 
point sources or from smoothly varying diffuse sources, 
such as intracluster and intrahalo dust. Thus the angular 
power spectrum should be a sensitive probe of the total 
dust content in the Universe. The existing estimes of 
the dust abundance from direct emission measurements 
make use of the sub-mm luminosity (e.g., Dunne, Eales 
& Edmunds 2003) or dust mass (e.g., Dunne et al. 2011) 
functions, they are generally based out of extrapolations 
of the measured bright galaxy counts. The anisotropy 
power spectrum should capture the integrated emission 
from faint sources, especially at the flux density scale 
that dominate the confusion noise. Separately, other es- 
timates of the cosmic dust abundance rely on the ex- 
tinction of optical light, especially with measurements 
that combine magnification and exti nction of quasars be- 
hind samples of foreground galaxies ([Menard et al.l[20Tol : 
[Ml) . It will be helpful to compare our direct emission 
measurement of the dust abundance with the extinction- 
based estimates since any differences can allow us to un- 
derstand the important of galaxies with hot dust that 
could be missed in SPIRE maps. We make use of a halo 
model to interpret the anisotropy power spectrum with 
the goal of measuring r2dust(-z), the cosmic abundance 
of dust relative to the critical density, as a function of 
redshift. 

To enable these measurements we make use of the wide 
field (~ 45 deg^) maps of Herschel-AThAS (Eales et al. 
2010) in the three GAMA areas along the equator, and 
select a single area that has the least Galactic cirrus 
confusion. This GAMA- 15 field involves 4 independent 
blocks of about 14 deg^, each overlapping with the adja- 
cent blocks by about 4 dcg^. We make use of the three 
overlap areas between the blocks to measure the power 
spectra at 250, 350 and 500 /im. The final power spec- 
trum is the average of the individual power spectra of 
each of the overlapping regions. While this forces us to 
make a measurement over a smaller area than the total 
survey area, our power spectrum measurement has the 
advantage that with two sets of cross-linked scans we 
can make independent measurements of the noise power 



spectrum. 

Our measurement approach is similar to that used for 
HerMES power spectra measurements (Amblard et al. 
2011; Viero et al. 2012) using multiple scans to gen- 
erate jack-knives of data to test the noise model. The 
total area used in HerMES measurements is about 12 
and 60 deg^, respectively, in Amblard et al. (2011) and 
Viero et al. (2012). However, H- ATLAS covers about 
120 deg^ in all three GAMA fields A measurement of the 
power spectrum in the whole of H- ATLAS GAMA areas 
requires an assumption about the noise power spectrum, 
since in regions with only one orthogonal scan or a sin- 
gle cross-linkeds scan, we are not able to separate the 
noise from the signal with data alone. In a future pa- 
per, we will present the power spectrum of the whole 
area using a noise model that is independently tested on 
various datasets to improve the confidence in separating 
noise in single cross-link scans. For now, wc make use 
of two cross-link scans for cross-correlations and auto- 
correlations to separate noise and sky signal. 

This paper is organized as follows. In Section [21 we 
briefly review how 250, 350 and 500 /im maps for the 
GAM A- 15 field were constructed using HIPE (Ott et all 
[2OT0h from raw time streams. In Section [3l we discuss 
how the auto and cross-correlation functions for each of 
the three fields were estimated, corrected, and assigned 
errors. The final power spectra are presented in Sec- 
tion [H The halo model used to fit the data and the 
luminosity function is discussed in Section [5l Finally, 
in Section [6| & [7| wc present our results and their im- 
plications, discuss future follow up work and give our 
concluding thoughts. 

2. MAP MAKING 

For this work we generate S PIRE maps using the 
MADmap (jCantalupo et al|[2009f ) algorithm that is avail- 
able within HIPE. The timeline data were reduced in- 
ternall y by the H-ATLA S team using HIPE version 
8.2.0 ( Pascale et al.l 1201 If ). The timelines were cali- 
brated with corrections applied for the temperature- 
drift and deglitched both manually and automatically. 
Astrometry corrections were also applied to the time- 
fines using offs ets between SDSS sources and the cross- 
identifications ([Smith et al.ll20lT[ ). In addition, a scan- 
by-scan baseline polynomial remover was applied to re- 
move gain variations leading to possible stripes. 

The map-maker, MADmap, converts the timeline data 
d{t) 

dit)=n{t) + Aip,t)xs{p), (1) 

with noise n{t) and sky signal s{p), given the pointing 
matrix A(p^ t) between pixel and time domain to a map 
by solving the equation 

m = {A^N-^A)-^A^N^^d . (2) 

Here N is the time noise covariance matrix and m is 
the pixel domain maximum likelihood estimate of the 
no iseless signal map given N and d. We refer the reader 
to lCantalupo et al ([2009( ) for more details of MADmap. 

The final maps we use for this work consist of four 
partially overlapping tiles, each containing two sets of 96 
scans in orthogonal directions (Fig. 1). The pixel-scale 
for the 250, 350 and 500 /im maps is 6", 8.333" and 12", 
respectively, corresponding to 1/3 of the beam size. 
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Fig. 1. — The Herschel- ATLAS GAMA-15 maps at 250 (top), 350 (middle) and 500 (bottom) fira with the three overlap regions used for 
the angular power spectrum measurements highlighted in dashed lines. 



In regions where the tiles do not overlap, the map at 
each wavelength consists of a single scan each in the two 
orthogonal scan directions. In the overlap region, we 
have two scans in each direction. As discussed below, we 
arc able to estimate the noise and signal power spectra 
independent of each other using the auto-correlations of 
the combined 4-scan map and the cross-correlations in- 
volving various jack-knife combinations. In Fig. 2 we 
show an example overlap region. 

3. POWER SPECTRA 



We now discuss the measurement of angular anisotropy 
power spectra in each of the three SPIRE bands. To 
be consistent with previo us measurements of t he SPIRE 
angular power spectrum (jAmblard et al.ll2011t) our maps 
are masked by taking a 50 mJy/beam flux cut and then 
convolving with the point response function. Such a flux 
cut, through a mask that removes the bright galaxies, 
also minimizes the bias coming from those bright sources 
by reducing shot-noise effects. The same mask also in- 
cludes a small number of pixels that do not contain any 
useful data, either due to scan strategy or data corrup- 
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Fig. 2. — The left overlap region in Fig. 1 at 250 Herschel-AThAS showing details of the background intensity variations without 
(top) and with (bottom) 5 > 50 mjy the bright source mask applied. This mask removes a substantial number of low-z bright galaxies 
detected in the areas used for the fluctuation study. 



tion. The combined mask removes roughly 13, 12 and 
15% of the pixels at 250, 350 and 500 /im, respectively. 
The fractions of masked pixels are substantially higher 
than the fractions of Amblard et al. (2011) of 1 to 2% as 
the ATLAS GAMA-15 field has a large density of z < 0.1 
spiral galaxies over its area relative to more typical ex- 
tragalactic fields used in the Amblard et al. (2011) study. 
These galaxies tend to be brighter, especially at 250 /im. 
While the fraction masked is larger, the total number of 
pixels used for this study is comparable to Amblard et 
al. (2011) with 2.9 X 10^, 1.5 x 10^ and 7.0 x 10^ at 250, 
350 and 500 /im in each of the three overlap regions. 

To measure the power spectrum in the final set of 
maps, we make use of 2D Fourier transforms. In gen- 
eral this is done with masked maps of the overlap re- 
gions, denoted Mi and M2 in real space. If we denote 
the 2D Fourier transform of each map as Mi and M2, the 
power spectrum, C;, formed for a specific / bin between 
between /-modes /i and I2, is the mean of the squared 
Fourier modes MiM| between li and 12- The same can 
be used to describe the auto power spectra, but with 
Ml = M2. 

The raw power spectra are summarized in Fig. 3. Here, 



we show the auto spectra in the total map, as well as the 
cross spectrum with maps made with half of the time- 
ordered data in each map. The difference of the two 
provides us with an estimate of the instrumental noise. 
At small angular scales (large £ values) the noise follows a 
white-noise power spectrum, with Ci equal to a constant. 
At large angular scales, the detectors show the expected 
l//-type of noise behavior, with the noise power spec- 
trum rising as Ci oc We fit a model of the form 



N,=N: 



1 



(3) 



and determine the knee-scale of the 1/f noise and the 
amplitude of noise power spectra. The noise values are 
No = 1.2 X 10^,5.3 X 10^ and 1.8 x 10^ Jy7sr at 250, 
350 and 500 ^m, comparable to the detector noise in the 
4-scan maps of the Lockman-hole used in Amblard et 
al. (2011). The knee at which 1/f noise becomes im- 
portant is Iq = 3730, 2920 and 3370, comparable to the 
expected knee at a wavenumber of 0.15 arcmin"^ given 
the scan rate and t he known properties of the detectors 
(|GrifRn et alllMol) . 
The raw spectra we have computed directly from the 
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Fig. 3. — The raw C; of the GAMA-15 overlap regions at 250 
(top), 350 (middle) and 500 (bottom) /im, respectively. The green 
points show the auto-power spectra computed from overlap regions 
using all 4 scans. This power spectrum is a combination of the real 
sky anisotropy power spectrum and the instrumental noise. We 
estimate the sky signal independent of noise by creating two sets 
of maps for each of the three overlap regions with two orthogonal 
scans each and then taking the cross power spectrum (red points) 
of those independent maps. The difference of these two spectra 
shows the instrumental noise power spectrum (black points). 



masked maps arc contaminated by several different ef- 
fects that must be corrected for. These issues are: the 
resolution damping from the instrumental beam, the fil- 
tering in the map-making process, and the fictitious cor- 
relations introduced by the bright source and corrupt 
pixel mask. Including these effects, we can write the 
measured power spectrum as 



(4) 



where C,' is the observed power spectrum from the 
masked map, B{1) is the beam function measured in a 
map, T{1) is the map-making transfer function, and Mu' 
is the mode coupling matrix resulting from the mask. 
Here, is the true sky power spectrum and is deter- 
mined by inverting the above equation. 

We now briefly discuss the ways in which we either 
determine or correct for the effects just outlined. 

3.1. The Map-Making Transfer Function 

Due to the finite number of detectors, the scan pattern, 
and the resulting analysis technique to convert timeline 
data into a map, the map we produce is not an exact 
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Fig. 4. — Map-making transfer function r(/) for the MADmap 
map making tool used for the GAMA-15 field anisotropy power 
spectrum measurement. The uncertainties in the transfer function 
are calculated from 100 random realizations of the sky as described 
in Section 3.1. 
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Fig. 5. — An example inverse-mode coupling matrix A-Z^j/ for one 
of the overlap regions (log scale) . 



representation of the sky. The modifications are associ- 
ated with the map-making process, relative to the true 
sky, are described by the transfer function T{1). We de- 
termine this by making 100 random realizations of the 
sky using Gaussian random fields derived from a first 
estimate of the H- ATLAS power spectrum. We sample 
those skies using the same timeline data as the actual ob- 
servations and analyze the simulated timelines with the 
same data reduction and map-making HIPE scripts for 
the actual data. We then compute the average of the 
ratio between the estimated power spectra and the input 
spectrum. This function is then the transfer function as- 
sociated with polynomial filtering and the map-making 
process. This transfer function, like the beam, represents 
a multiplicative correction to the data. We divide the 
estimated power spectrum of the data by this transfer 
function to remove the map-making pipeline processing 
effects. 

In Fig. 4 we show the transfer functions at 250, 350 
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and 500 with 68% error bars taken from the stan- 
dard deviation of 100 simulations. The transfer function 
is such that it turns over from 1 at both large angular 
scales, corresponding to roughly the scale of an individ- 
ual scan length, and the beam scale. The large-scale 
deviation, which is wavelength independent, is due to 
the polynomial removal from each timeline of data, while 
the turnover at small angular scales is due to the cut-off 
imposed by the instrumental point response function or 
the beam. The transfer function is more uncertain at the 
large angular scales due to the finite number of simula- 
tions and the associated cosmic variance resulting from 
the field size. Given its multiplicative nature, errors from 
this transfer function are added in quadrature with rest 
of the errors. 

3.2. The Beam 

Following Amblard et al. (2011), the beam function is 
derived from Neptune observations of SPIRE. The Nep- 
tune timeline data are analyzed with the same pipeline 
and our default map maker in HIPE. The resulting beam 
functions are similar to those of Amblard et al. (2011) 
and we find no detectable changes resulting from the two 
different map rnakers between this work and the SMAP 
(jLevenson et al.ir2010D pipeline of the SPIRE Instrument 
Team used in Amblard et al. (2011) and Viero et al. 
(2012). This is primarily due to the fact that the beam 
measurements involve a large number of scans and Nep- 
tune is several orders of magnitude brighter than the 
extragalactic confusion noise. Wc interpolate the beam 
function measured from Neptune maps in the same i 
modes at which we compute our anisotropy power spec- 
tra. This beam transfer function, B{1), at each of the 
wavelengths represents a multiplicative correction to the 
data. Similar to Amblard et al. (2011), we compute the 
uncertainty in the beam function by computing the stan- 
dard deviation of several different estimates of the beam 
function by subdividing the scan data to 4 different sets. 
The error on the beam function in Fourier space is prop- 
agated to the final error and is added in quadrature with 
rest of the errors. 

3.3. Mode-Coupling Matrix 

The third correction we must make to the raw power 
spectrum involves the removing of fictitious correlations 
between modes introduced by the bright sources and con- 
taminated or zero-data pixel mask. Due to this mask the 
2D Fourier transforms are measured in maps with holes 
in them. In the power spectrum these holes result in a 
Fourier mode coupling that biases the power spectrum 
lower at large angular scales and higher at smaller an- 
gular scales. This can be understood since the modes at 
the largest angular scales, like the mean of the map, are 
broken up into smaller scale modes with any non-trivial 
mask. 

To c orrect for the mask we make use of the method 
used in iCoorav et al.l (|2012| ). The method involves cap- 
turing the effects of the mask on the power spectrum into 
a mode-coupling matrix Mw . The inverse of the mode 
coupling matrix then removes the contamination and cor- 
rects the raw power spectrum to a power spectrum that 
should be measurable in an unmasked sky. The correc- 
tion both restores the power back to the large angular 



scale modes by shifting the power away from the small 
angular scale modes, especially those at the modulation 
scale introduced by the mask. 

To generate Mw we apply the mask to a map consisting 
of a Gaussian realization of a single Z-modc and take the 
power spectrum of the resulting map. This power spec- 
trum represents the shuffling of power the mask performs 
on this specific /-mode among the other Z-modes. This 
process is repeated for all /-modes and these effects of 
the mask on each mode are then stored in a matrix. This 
matrix, Mw now represents the transformation from an 
unmasked to a masked sky by construction. By inverting 
this matrix, shown in Fig. 5, we arc left with the trans- 
formation from a masked to an unmasked sky removing 
the fictitious couplings induced by the mask applied to 
the raw power spectra. The matrix Mw behaves such 
that in the limit of no /—mode coupling Mw = fsky^w 
where /sky is the fraction of the sky covered. Thus in 
the limit of partial sky coverage the correction becomes 
the standard formula with C/ = /skyC';- For more de- 
tails, including figures demonstrat hig the robustness o f 
the method, we refer the reader to ICoorav et al.l ()2012[ ). 

4. POWER SPEGTRUM RESULTS 

The final power spectrum Ci at each of the three wave- 
lengths is shown in Fig. 6 (left panels). The final error 
bars account for the uncertainties associated with the (a) 
beam, (b) map making transfer function, (c) instrumen- 
tal or detector noise (Fig. 3), and the cosmic variance 
associated with the finite sky coverage of the field. In 
Fig. 6 we compare these final H- ATLAS GAMA-15 power 
spectra with measurements of the GFIRB anisotropy 
power spectrum measurements from HcrMES (Amblard 
et al. 2011) and Planck (Planck collaboration 2011) 
team measurements. We find general agreement, but 
we also find some differences. At 250 /im, wc find the 
amplitude to be larger than the existing SPIRE mea- 
surements of the power spectrum at 250 fj,m in HerMES 
while the measurements are more consistent at 350 and 
500 iim. We attribute this increase to the wide coverage 
of H-ATLAS and the presence of a large surface den- 
sity of galaxies at low redshifts. While most of these 
galaxies are masked we find that the fainter population 
likely remains unmasked and contributes to the increase 
in the power that we have seen. This conclusion is also 
consistent with the strong cross-correlation between de- 
tected SPIRE sources in GAMA fields of H-ATLAS and 
the SDSS redshift survey (e.g., Guo et al. 2011). The 
difference between Herschel-SPIRE measurements and 
Planck measurements are discussed in Planck collabora- 
tion (2011) and we refer the reader to that discussion. 
We continue to find differences between our measure- 
ments and Planck power spectra at 350 yum, even with 
Planck data corrected for the frequency differences and 
other corrections associated with the source mask, as dis- 
cussed in Planck collaboration (2011). 

Note that the power spectra in the left panels of Fig. 6 
asymptote to a C; constant. This is the shot-noise 
coming from the Poisson behavior of the sources. In 
Fig. 6 (right panels) we show the final power spectra 
plotted as 1'^Ci/2tt, with the Poisson noise removed at 
each band. They now reveal the underlying clustering of 
submillimeter galaxies. With sources masked down to 50 
mJy, our shot-noise amplitudes are 6700±140, 4400±130 
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Fig. 6. — Final angular power spectra of CFIRB anisotropies in the H-ATLAS GAMA-15 field at 250 (top), 350 (middle) and 500 (bottom) 
fira. In the left panels, the power spectra are plotted as C; prior to the removal of the shot-noise term. Here we compare the power spectra 
measured with H-ATLAS data to Planck and previous Herschel results from HerMES. In the right panels wc show the power spectra as 
PCi/2tt after removing the shot-noise level at each of the wave bands. We add the uncertainty associated with the shot-noise level back 
to the total error budget in quadrature. This results in the increase in errors at high multipoles or small angular scales. The curves show 
the best-fit model separated into 1 and 2-halo terms (see text for details) and the total (orange line). The solid line that scales roughly 
as l^Ci ~ is the best-fit Galactic cirrus fiuctuation power spectrum. Due to the high cirrus fiuctuation amplitude and clustering, 

the H-ATLAS power spectrum in the GAMA-15 field at 250 fiia is higher than the existing HerMES results, while the measurements are 
generally consistent at 350 and 500 /im. 



and 1900±90 JyVsr at 250, 350 and 500 fim, respectively 
(see Table 1). We determine the Poisson noise uncertain- 
ties based on the overall fit to C; measurements at the 
three highest £-bins. 

For comparison to our shot-noise values, the shot-noise 
values of Amblard et al. (2011) are 6100± 120, 4600±70 
and 1800 ± 80 JyVsr at 250, 350 and 500 /im, respec- 
tively. While the shot-noise values are consistent at 350 
and 500 //m, we find an increased shot-noise amplitude at 
250 /im, consistent with the higher amplitude of the clus- 
tering part of the power spectrum. In addition to Planck 
and Amblard et al. (2011) HerMES measurements, in 
Fig. 6 (right panels) wc also compare our measurements 
to more recent Viero ct al. (2012) HerMES measure- 
ments. At 350 and 500 /im the difference between all of 
Herschel-SPIRE measurements and Planck is clear. 

At 250 /xm, we find that our measurements have a 



higher amplitude at all angular scales relative to previous 
SPIRE measurements. At large angular scales, we find 
that the increase is coming from the higher intensity of 
cirrus in our GAMA-15 fields (Bracco et al. 2011). The 
cirrus properties as measured from the power spectra are 
discussed in Section 6.1. As part of the discussion related 
to our results on the galaxy distribution that is contribut- 
ing the far-IR background power spectrum (Section 6.2), 
wc will explain the difference between the HerMES and 
H-ATLAS power spectrum at 250 /im as due to an excess 
of low-rcdshift galaxies in the H-ATLAS GAMA-15 field 
(Rigby et al. in prep). The measurements shown in Fig. 6 
right panels constitute our final CFIRB power spectrum 
measurements in the H-ATLAS GAMA-15 field. These 
power spectra values are tabulated in Table 2. We now 
discuss the model used for the interpretation leading to 
the best-fit model lines shown in Fig. 6 (right panels). 
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5. HALO MODELING OF THE CFIRB POWER 
SPECTRUM 

To analyze the H- ATLAS GAMA-15 power spectra 
measurements we implement t he conditional luminos- 
ity fu n ction (CLF) approach of [Giavalisco fc DickinsonI 
((20MI) : lLee"enil (12001 i fPe Bernardis fc Cooravl (12011) 7 
We recall below the main features of the model and re- 
fer the reader to these works for more details. The goal 
is to work out the relation between IR luminosity and 
halo masses of the galaxies that are contributing to the 
CFIRB power spectrum. We populate halos with the 
best-fit Lir(M) relation from the data and use that to 
determine the abundance of dust (Jldust) in the Universe. 
The CLF approach proposed here improves over several 
assumptions that were made in Amblard et al. (2011) 
to interpret the first Herschel-SPIRE anisotropy power 
spectrum measurements. 

First, the probability density for a halo or a sub-halo of 
mass M to host a galaxy with IR luminosity L is modeled 
as a normal distribution with: 



P{L\M) 



1 



2'n(JL{M) 



cxp 



{L-L{M)f 



(5) 



The relation between the halo mass and the average lu- 
minosity L{M) is expected to be an increasing function 
of the mass w ith a characterist ic mass scale Mqi and we 
can write (^see lLee et al.l ()2009[ )) 



m)^L.[—) exp 



M 



-A 



(6) 



As already discussed by iLee et ahl (|2009[ ) these param- 
eterizations do not have a specific physical motivation, 
except for the requirement that the luminosity increases 
as an increasing function of the halo mass and offer the 
advantage that one can explore a large range of pos- 
sible shapes for the luminosity-mass relation. While 
there is no motivation to use this specific form over an- 
other, certain models of galaxy formation do predict a 
L{M, z) relation and our results based out of the model- 
fits to CFIRB power spectrum can be compared to those 
model predictions. In particular, the model of Lapi et 
al. (2011) predicts L{M,z) oc M(l + z)^-^, while the 
cold-flow accretion model of Dekel et al. (2009) predicts 
L(Af,z)cxMi-i5(l + z)2-25. 

The total halo mass function is given by the number 
density of halos or sub-halos of mass M. The contribu- 
tion of h alos nh(M) is taken to b e the Sheth & Tormen 
relation (jSheth fc Tormeiil Il999( ) . The sub-halos term 
can be modeled through the number of sub-halos of mass 
m inside a parent halo of mass Alp, N{m\Mp). The total 
mass function is then written as 



nriM) ^ nh{M) + nsh{M) , 
where nsh{M) is the sub- halo mass function 

nsh{M)^ I N{M\Mp)nh{Mp)dMp. 



(7) 



(8) 



Here we paramet erize N(m\M) following the semi- 
analytical model of Ivan de Bosch et al.l ()2005f ). 

Neither the normalization nor the slope of the sub-halo 
mass function are universal and both depend on the ratio 



TABLE 1 

Parameter values from MCMC fits to the H-ATLAS 
gama-1.5 angular power spectra at 250, 350 and 500 ^m. 



HOD 




n fiq -1- n 04 




Bi 

HI 


0.09 ± 0.05 




loP'f Tyti i T 

logrAfn/jV/o) 

Pm 


9.52 ± 0.08 
11.5 ± 1.7 
-2.9 ±0.4 


CFIRB SED 


Tdust 


37±2 K 




/3dust 


unconstrained 


Cirrus 


^-^250 


3.5 ± 1.3 X 10^ Jy^/sr 




'-'350 

r'i=230 

'-'500 

cirrus 


1.2 ± 1.0 X 10* Jy^/sr 
1.1 ± 0.9 X 10^ Jy2/sr 
21.1 ± 1.9 K 




/^cirrus 


2.9 ±0.8 


Poisson 


SAf250 


6700 ± 140 Jy^/sr 




5Af350 


4400 ± 130 Jy2/sr 




SNrMO 


1900 ±90 Jy2/sr 



between the parent halo mass and the non-linear mass 
scale, M*. is defined as the mass scale where the 
rms of the density field tT(M, z) is equal to the critical 
over-density required for spherical collapse 5c{z). The 
contribution of central galaxies to the halo occupation 
distribution (HOD) is simply the integral of P{L\M) over 
all luminosities above a certain threshold Lq, either fixed 
by the survey or a priori selected so that 



(^c(A/)) 



L>L„ 



P{L\M)dL, 



(9) 



which, in the absence of scatter, reduces to a step func- 
tion Q{M—Mq), as expected. Note that all integrals over 
the luminosity L also have a redshift-dependent cut-off 
at the upper limit corresponding to the flux cut of 50 
mjy that we used for the power spectrum measurement. 
For the satellite galaxies, the HOD is related to the sub- 
halos 



L>L„ 



dL / dmN{m\M)P{L\m) 



(10) 



The total HOD is then 

(Mot(M))L>L„„„ = {Nh{M))L>L„ 



{Nsh{M)) 



L>L„ 



(11) 



We account for the possible redshift evolution of the 
luminosity-halo mass relation by introducing the param- 
eter Pm and rewriting the mass scale A/q; as: 



(12) 



Under the assumption that the central galaxy is at the 
center of the halo and that the halo radial profile of satel- 
lite galaxies within dark matter halos follow that of the 
dark matter, gi ven by the Navarro , Frenk and White 
(NFW) profile (jNavarro et al.|[T997D . we can write the 
1-halo and 2-halo terms of the three-dimensional power 
spectrum: 

P^'^ik) = 4 / dM{NT{NT - l))u{k,MYnH{M){li) 
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Fig. 7. — Best-fit halo occupation distribution and the Icr range 
at 2 = 1 for three cases involving Lin > 10'',10^'' and 10^^ Lq. 
The three lines to the top show the different power-laws for compar- 
ison with the shape of the HOD. The satellite galaxies contribution 
has a slope ~ 1 when Ljj)^ ~ 10^ Lq. 



Fig. 9. — Normalized redshift distributions of FIR-bright galax- 
ies predicted by our model for two different flux density cuts at 250 
^l (thick solid line for 1 mjy < 5 < 10 mjy and thick dashed line 
for 10 mjy < S < 50 mjy). For comparison in corresponding thin 
linos wo show the measured redshift distributions for the Herschel- 
selected g alaxies (HSG s ) at th e same flux density bins with optical 
spectra in lCasev et al.l l|2012l '). 
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Fig. 8. — L uminosity fu n ctions predicted by our mode l compared 
to dat a from lEales et all POIOD (0.2 < z < 0.4) and ILapi et al.1 
H20111 'l (1.2 <z < 1.6,2 < 2 < 2.4). The shaded region corresponds 
to the 68% confidence level. 



where u{k,AI) is the NFW profile in Fourier space and 
Hg is the galaxy number density 



= / dM{Ng{M))nh{M). 



(14) 



The second moment of the HOD that appear in eq. [13] 
can be simplified as 



{Nt{Nt - 1)) ^ {Nrf - (Nh)' 



(15) 



and the power index p for the NFW pro file is p = 1 whe n 
{NriNr - 1)) < 1 and p = 2 otherwise (|Lee et al.ll2009l ). 
The two-halo term of galaxy power spectrum is 



P2''(fc): 



— / dM{NT{M))u{k,M)nh{M)b{M) 



xPiin(fc) , 



(16) 



where Punik) is the linear power spect rum and b(M) is 
the linear bias factor calculated as in (|Coorav fc ShethI 
|2002| ). The total galaxy power spectrum is then Pg(fc) = 
P^f'{k) + P^^{k). 

The observed angular power spectrum can be related 
to the three-dimensional galaxy power spectr um through 
a redshift integration along the line of sight (jKnox et alJ 
[200lHA^ard fc Cooravll2007l) : 



dz 



dx 
dz 



U^)j,,iz)Pgi£/x,z),{l7) 



where x isjthe comoving radial distance, a is the scale 
factor and ji,{z) is the mean emissivity at the frequency 
1/ and redshift z per comoving unit volume that can be 
obtained as: 



dL(f){L, z)- — , 
47r 



(18) 



Here the luminosity function is 



0(i, z)dL ^dL j dMP{L\M)nT{M, z) . (19) 

To fit data at different frequencies we assume that the 
luminosity-mass relation in the IR follows the spectral 
energy distribution (SED) of a modified black-body (here 
we normalize at 250fim at z = 0) with 

L,(M) = i2.5o(M)(l - e-^)B{i^o,Ta) , (20) 

where Td is the dust temperature and the optical depth 
is T = {^)^'' when B{vQ^Td) is the Planck function and 
L25o{M) is given by eq dH). 

The final power spectrum is a combination of galaxy 
clustering, shot-noise and the Galactic cirrus such that 

^tot ^ c-pFIRB ^ (^cirrus ^ (jSn ^ ^^^^^ ^CFIRB 

power spectrum derived above and Cf^ is the scale- 
independent shot-noise. To account for the Galactic cir- 
rus contribution to the CFIRB, we add to the predicted 
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Fig. 10. — Best-fit determination of mean emissivity at 250 /im as 
a function of the redshift, uj^{z) (thick soUd line), and its la error 
from the MCMC model fits (grey shaded region) for sources with 
^250 < 50mjy. We show several model predictio ns from the litera- 
ture HValiante et aLllSOOglTBethermin etaIl f20TTh and compare our 
estimates to the determinations from the halo mode l fits to the 
CFIR B power spectra by [Smblard et al. (2011,) and IViero et al.l 
II2012I '). lAmbiard et al.l 11201 ll ) measurements involve a binned de- 
scription of j^jz) with la errors determined from the fit, while 
IViero et al.l lj2012l ') result is the best-fit relation for their work. 



angular power spectrum a cirrus power-l aw power spec- 
trum w ith the same shape of that used bv lAmblard et aP 
(j2011|) . where the authors assumed the same cirrus 
power- law power-s pectrum shape from measurements of 
IRAS and MIPS (|Lagache et alj I2007D at 100 ^im with 
Ci (X r" with n = 2.89 ± 0.22. In Amblard et al. 
(2011) this 100 /xm spectrum was extended t o longer 
wavele ngths using the spectral dependence of ISchlegell 
(|T99l . Here we rescale the amplitude of the cirrus power 
spectrum with amplitudes C""^^ at each of the three 
wavelengths {i = 250, 350 and 500 iJ.) taken to be free pa- 
rameters and model-fit those three parameters describing 
the amplitude as part of the global halo model fits with 
the MCMC approach. 

6. RESULTS AND DISCUSSION 

We fit the halo-model described above to the 250, 350 
and 500 fim CFIRB angular power spectrum data for the 
H- ATLAS GAMA-15 field by varying the halo model pa- 
rameters and the SED parameters. The dimension of the 
parameter space is thus 12 with free parameters involv- 
ing Td, Pd, ai, (3i, Lq, pm, C^'^'i-^so and SNi. We 
make use of a Markov Chain Monte Carlo (MCMC) pro- 
cedure, modified from the publicly available CosmoMC 
([Lewis fc Bridld [2002D ■ with a conver gence diagnostics 
based on the Gelman- Rubin criterion (IGelman &: Rubiiil 
I1992D . To keep the number of free parameters in the 



halo model manageable, we a priori constrain the Mg; in 
equation [12] to the value of log Mqi / Mq = 11. 5± 1.7 as 
determined by a fit to the low- redshift lumin osity func- 
tion at 250 fim ([De Bernardis &: Cooravl2012l) using data 
from Vaccari et al. (2010) and Dye et al. (2010). The 
best fit parameters and the uncertainties from the halo 
model fits are listed in Table [TJ 

6.1. Cirrus Amplitude and Cirrus Dust Temperature 




36 36.5 

T,(K) 

Fig. 11. — 68% and 95% confidence level constraints on T^j and 

/3dust- 



We now discuss some of the results starting from our 
constraints on the cirrus fiuctuations. The cirrus ampli- 
tudes have values of (3.5 ± 1.3) x 10^ (1.2 ± 1.0) x lO" 
and (1.1 ± 0.9) x 10^ Jy^/sr at 250, 350 and 500 urn, 
respectively, at € = 230 corresponding to 100 arcminute 
angular scales. These values are comparable to the cirrus 
amplitudes in the Lockman-Hole determined by Amblard 
et al. (2011). The GAMA-15 area we have used for this 
study is thus comparabale to some of the least Galactic 
cirrus contaminated fields on the sky. For comparison, 
the GAMA 9 hour area studied by Bracco et al. (2011) 
has cirrus amplitudes of ~ 3 x 10'', 2 x 10^ and 1 x 10^ 
at 250, 350 and 500 /x, respectively. These are roughly 
a factor of 100 larger than the cirrus fluctuation ampli- 
tude in the GAMA-15 areas used here. The third field 
we considered for this study in GAMA 12 hour area was 
found to have cirrus amplitudes that are roughly a factor 
of 20 to 30 larger. 

In order to determine if the cirrus dust in the GAMA- 
15 field is comparable to dust in the high cirrus intensity 
regions such as the GAMA-9 field, we fitted a modified 
blackbody model to the cirrus rms fluctuation amplitude. 
We found the dust temperature and the dust emissivity 
parameter /? to be 21.1±1.9 K and 2.9±0.8, respectively. 
The results from the same analysis at 100 arcminute-scale 
rms fluctuations are 20.1 ±0.9 and 1.3±0.2 for dust tem- 
perature and emissivity, respectively. Even though the 
cirrus amplitude is lower with rms fluctuations, yXSf^^, 
at a factor of 10 below the GAMA-9 area studied in 
Bracco et al. (2011), we find the dust temperature to 
be comparable. It is unclear if the difference in the dust 
emissivity parameter is significant or captures any phys- 
ical variations in the dust from high to low cirrus inten- 
sity, especially given the well-known degeneracy between 
dust temperature and /3. Fluctuation measurements in 
all of the 600 sq. degree H- ATLAS fields should allow a 
measurement of /3 as a function of cirrus amplitude. 

6.2. Faint Star-Forming Calaxy Statistics 

Moving to the galaxy distribution, in Fig. 7, we show 
the halo occupation distribution at z = 1 correspond- 
ing to the best-fit values of the parameters and the Icr 
uncertainty region for three different luminosity cut-off 
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Fig. 12. — The cosmic density of dust O^uat vs redshift as determined from the CFIRB power spectra from H-ATLAS GAMA-15 field 
(shaded region). The thickness of the region corresponds to the 1 — cr ranges of the halo model parameter uncertainties as determined by 
MCMC fits to the data (Table 1). We also compare our estimate to previous measurements in the literature. The measurements labeled 
H-ATLAS dust mass function are from the low-rodshift dust mass function measurements in Dunne et al. (2011). The other estimates are 
based on extinction measurements from the SDSS (e.g., Menard et al. 2010, 2012; Fukugita 2011; Fukugita & Peebles 2004) and 2dF (e.g.. 
Driver et al. 2007). 



values. At z = 1, as show in Fig. 7, for ijR > 10^ 
Lq galaxies, the HOD drops quickly for masses smaller 
than \og{Mjnin/MQ) ~ 10.7 and the high-mass end has 
a power- law behavior with a slope ^ 1. By design, this 
halo model based on CLFs has the advantage that it does 
not lead to unphysical situations with power-law slopes 
for the HOD greater than one as found by Amblard et 
aL (2011). 

Both the HOD and the und erlying luminosity-mass 
relatio ns are consistent with iDe Bernardis fc Cooravl 
(j2012| ). where a similar model was used to reinterpret 
Amblard et al. (2011) anisotrop y measurement. The key 
differe nce between the work of IDe Bernardis fc Cooravl 
(|2012l) and the work here is that we introduce a dust 
SED to model-fit power spectra measurements in the 
three wavebands of SPIRE, while in earlier work only 
250 /im measurements were used for the model fit. For 
comparison with recent model descriptions of the CFIRB 
power spectrum, we also calculate the effective halo mass 
scale given by 



Afeff = / dMnh{M)M 



Nt{M) 



(21) 



Msiin a t z 



With this definition we find M^s = 3.2 x 10^ 

2, con siste nt with the effective mass scales of lShang et al.l 
(|20Tl[) and IDe Bernardis fc CooTavl (|20T2| ) of M^s ~ 4 x 
10^^ and slightly lower than the value of ~ 5 x 10^^ from 



IXia et all (l20Tl) . 

The MCMC fits to the CFIRB power spectrum data 
show that the charicteristic mass scale Afpz evolves with 
redshift as (1 -I- z)"^-^^" "'. In order to compare this 
with existing models, we convert this evolution in the 
charachteristic mass scale to an evolution of the i(M, z) 
relation. As L{M) cx (M/Mq;)""', we find L{M,z) oc 
M;"(l + z)-P"°". Using the best-fit values, we find 
L(M,z) (X MO ™±° 05(1-Hz)2-0±o.4_ inLapietaL (2011), 
their equation 9 with the SFR as a measure of the IR lu- 
minosity, this relation is expected to be M(l -I- z)^-^. In 
Dekel et al. (2009), the expectation is M'^ '^^{\ + zf-"^^. 
While we find a lower value for the power-law depen- 
dence on the halo mass with IR luminosity, the redshift 
evolution is consistent with both these models. 

To test the overall consistency of our model relative to 
existing observations at the bright-end, in Fig. 8, we com- 
pare the predicted luminosity functions 250 /im-selected 
galaxies in several redshift bins with existing measure- 
ments in the literature from Eales et al. (2010) and Lapi 
et al. (2011). The former relics on the spectroscopic 
redshifts in GOODS fields while the latter makes use of 
photometric redshifts. We find the overall agreement to 
be adequate given the uncertainties in the angular power 
spectrum and the resulting parameter uncertainties of 
the halo model. In future, the overall modeling could 
be improved with a joint-fit to both the angular power 
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spectra and the measured luminosity functions. 

In Fig. 9 we show the predicted redshift distributions 
of the 250 /im-selected galaxies in two 250-/im flux den- 
sity bins in our model with a comparison to a measured 
redshift distribution with close to 900 optical spectra of 
Herschel-selected galaxies with Keck/LRIS and DEIMOS 
in Casey et al. (2012). While there is an overall agree- 
ment, for the brighter flux density bin, the measured 
redshift distribution shows a distinct tail a small, but 
non-negligible, fraction of galaxies at z > 2. It is unclear 
if those redshifts suggest the presence of bright galaxies 
that are lacking in our halo model or if tho se redshifts are 
associated with lensed sub -mm galaxies (jNegrello et al.l 
l2010t IWardlow et al.|[2012[) with intrinstic fluxes that are 
below 10 mJy. If lensed, due to magniflcation boost, such 
fainter galaxies will appear in the brighter bin. We also 
note that the current halo model ignores any lensing ef- 
fect in the anisotropy power spectrum. Existing models 
suggest that the lensing rate at 250 ^m with flux den- 
sities below 50 mJy is small. At 500 ^m, however, the 
lensed counts are at the level of 10% (Wardlow et al. 
2012). While we do not have the signal-to-noise ratio for 
a lensing analysis of the far-IR background anisotropics 
with the current data and the power spectrum, a future 
goal of sub-mm anisotropy studies must involve charac- 
terizing the lensing modification to the power spectrum. 

In Fig. 10 we show the redshift evolution of the emis- 
sivity predicted by our model at 250 /im according to 
equation 1181 The shaded region shows the la uncer- 
tainty associated with the best-fit model. For comparison 
we show the results oflViero et al.l (l2012D;lValiante et al.l 
(l2009t): iBethermin et al.l (|2011[ ): lAmblard et al.l (|20lH ): 
iGispert et al.l ( 2000[ ). The distribution predicted by our 
fit is consistent for a wide ran ge of redshifts (up t o 
z > 3) with iViero et all (l2012fl: iVahante et all (|2009( ): 
IBethermin et all (120111^ " lAmblard et al.l (|2011| ). The re- 
cent fit of IViero et all (|2012f ) to the HerMES angular 
power spectra shows a lower emissivity at both low- 
redshift {z < 0.5) and high-redshift {z > 2.5) ends. 

The excess in the emissivity at the low- redshifts (z < 
0.1) partly explains the difference in the power spectrum 
amplitude at 250 fim between the previous angular power 
spectra and H-ATLAS data. As discussed earlier, the 
GAMA-15 field of H-ATLAS is known to contain an over- 
density of low-redshift galaxies. The brightest of these 
sources with 5*250 > 50 mJy are clearly visible in Fig. 2 
when comparing the original and masked maps. While 
the 5250 > 50 mJy mask is expected to remove a sub- 
stantial fraction of the low-z population, we expect a 
fraction of the fainter ones to remain. Such galaxies are 
not present in the well-known extragalactic fields of the 
HerMES survey such as Lockman-Hole and the NDWFS- 
Bootcs field. The difference is a factor of ^ 2 amplitude 
increase in the power spectrum at 250 ^m. As the excess 
population is primarily at low redshifts, the difference 
only shows up at 250 ^m, while we do not see any sig- 
nificant difference at 350 and 500 fj,m between HerMES 
and H-ATLAS power spectra. The final result of this is 
to increase the emissivity at 250 ^m at lowest redshifts 
z < 0.1 in our model r elativ e to the emissivity function 
derived in IViero et all ()2012D . 

The H-ATLAS GAMA-15 field also shows an overall in- 
crease of bright counts at 250 /im relative to the HerMES 
fields (Rigby et al. in prep) and we verified that our sug- 



gestion of a factor of 2 increase in the power spectrum 
is coming from low-redshift galaxies is consistent with 
the differences in the number counts. The difference in 
the counts also explains the increase in the shot-noise 
at 250 fim relative to the value found in HerMES power 
spectra. These differences generally suggest that large 
field-to-field variations in the angular power spectrum 
with variations well above the typical Gaussian cosmic 
variance calculations. Such variations are readily vis- 
ible when comparing individual field power spectra in 
IViero et al. , (2012) (their Fig. 3). 

Through our joint model-fit to 250, 350 and 500 yum 
power spectra, we also determine the SED of far-IR back- 
ground anisotropics. To keep the number of free parame- 
ters in our model small, here we assume that the SED can 
be described by an isothermal black-body model. The 
best-fit dust temperature value that describes the far-IR 
fluctuations is 37 ± 2 K while the emissivity parame- 
ter /3 is unconstrained. In Fig. 11 we show the best- fit 
68% and 95% confidence level intrvals of Tj, and /3, af- 
ter marginalizing over all other parameters of the halo 
model. This figure makes it clear why we are not able to 
determine /3 with the current data due to degeneracies 
between the model parameters. The dust temperature 
we measure should be considered as the average dust 
temperature of all galaxies that is contributing to far-IR 
background anisotropy power spectrum. The dust tem- 
perature is higher than the typical 20K dust temperature 
derived from the absolute background spectra at far-IR 
wavelengths from exp eriments such as FIRAS and Planck 
(|Lagache et al.ll2000l) . 

This difference in the dust temperature could be un- 
derstood since the absolute measurements, especially at 
degree angular-scale beams, are likely to be dominated 
by the Galactic cirrus, and thus the temperature mea- 
surement could be biased low. The dust temperature 
wc measure from the far-IR power spectra is fully con- 
sistent with the the value of 44 ± 7K by Shang et al. 
(2011) in their modeling of the Planck far-IR power spec- 
tra (assuming the fixed value /3 = 2). Thus, while the 
best-fit SED model of the absolute GIB may suggest a 
low temperature value, the anisotropics from approxi- 
mately 1-30 arcminute angular scales follow a SED with 
a higher dust temperature value. Separately, we also 
note that our dust temperature of 37 ±2 K is also con- 
sist ent with the average dus t temperature valu es of 36±7 
K (jChapman et al. ,200 5; Dunne et al.l 120001 ) for high- 
z SCUBA-selected sub-mm galaxies, but is somewhat 
higher than the average dust temperatur e value of 28 zt 8 
K for Herschel-selected bright galaxies in lAmblard et al] 
(|2010f) . The Amblard et al. (2010) value is dominated 
by low-redshift (z ^ 0.1) galaxies with Herschel identi- 
fications to SDSS redshifts. In the local Universe, most 
dusty late-type gala xies show cold dust with tempera- 
tures around 20 K (|Galametz et al.ll2012t iDavies et al.l 
120121 ). The higher temperature we find for the far-IR 
background anisotropics then suggests that the average 
interstellar radiation field in galaxies at z ~ 1 to 2 that 
dominate the dust emissivity is higher by a factor of 2^ 
when compared that local late-type galaxies. 



6.3. Cosmic Dust Abundance 
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The model described above allows us to estimate the 
fractional cosmic dust density: 



f^dust(^) 



1 



dLcj){L,z)Mdust{L), (22) 



where Mdust is the dust mass for a given IR luminosity 
and po is the critical density of the Universe. Here we 
make use (/)(L, z) as derived by the halo model fits to the 
far-IR background power spectra. 

To conv ert luminositie s to dust mass, we follow equa- 
tion 4 in (iFu et al.ll2012f ). This requires an assumption 
related to the dust mass absorption coefficient, Kd- It is 
generally assumed that the opacity follows (x 
with a normalization of Kd = 0-07 ib . 02 ve? kg~^ at 850 
^m (jPunne et all 120001 : I James et al1l2002D . This nor- 
malization, unfortunately, is highly uncertain and could 
easily vary by a factor of few or more (see discussion 
in James et al. 2002). The value we adopt here is ap- 
propriate for dusty galaxies and matches well with the 
integrated spectrum of the Milky Way. 

The conversion to dust mass also requires the SED of 
dust emission. Here we make use of the average dust tem- 
perature value of 37±2 K as determined by the model-fits 
to the angular power spectra. As (3 is undetermined from 
the data, we take its range with a prior between 1 and 2.5, 
consistent with typical values of 1.5 or 2 that is generally 
assumed in the literature. When calculating f2dust(^) we 
marginalize over all parameter uncertainties so that we 
fully capture the full likelihood from the MCMC chains 
given the prior on /3. Note that our assumption of a con- 
stant dust temperature is at odds with local late-type 
spirals that show much lower temperatures. However, 
it is also known that there are some sub-mm galaxies, 
especially those that are radio bright, with dust temper- 
atures in the excess of 60 K. Thus with a value of 37 ±2 
K we may be using a representative average value for the 
dust temperature and an average of the dust SED for all 
galaxies at a variety of redshifts. Finally there are some 
indications that the dust temperature is IR luminosity 
dependent (see the discussion in Amblard et al. 2010). 
If that remains to be the case then the correct approch 
with eq. [21] will be to take into account that luminosity 
dependence as seen in the observations. Given that the 
current indications are coming from small galaxy sam- 
ples, we do not pursue such a correction, but highlight 
that future studies could improve our dust abundance 
estimate. 

In Fig. 12 we show our results. In addition to the di- 
rect emission estimate that we have considered here, we 
also show the low-z dust abundances by integrating over 
the dust mass functions in Dunne et al. (2011). Those 
dust mass functions are limited to 2 < 0.5 due to the lim- 
ited availability of spectroscopic data at higher redshifts. 
While in principle dust mass function captures the to- 
tal dust of detected galaxies, the mass functions can be 
extrapolated to the faint-end, as has been done here, to 
acocunt for the fainter populations below individual de- 
tection levels. Thus the abundances from mass functions 
must agree with the estimates based on the anisiotropy 
measurements. We do not use our halo model to estimate 
the dust abundance at z < 0.05 since our halo model is 
normalized to the luminosity function of dusty galaxies 
at low redshifts. 



Our measurement indicates that the dust density 
ranges between fidust — 10~^ to 8 x 10~^ in the red- 
shift range z = 0.5 — 3. We note that the fidust pre- 
di ction of this work has a smaller uncertainty than that 
in IDe Bernardis fc Cooravl l|2012D where the estimation 
was done assuming a larger range for Td and Pd- In 
equation mi we integrate over luminosities L > 10^ Lq. 
However in this calculation the choice of minimum lumi- 
nosity is less relevant, since the uncertainty on the dust 
density estimate is dominated by the large uncertainties 
of dust temperature and spectral emissivity index f3. 

Fig. [T^ also summarizes the dust-density rneasure - 
ments oflFukugita fc Peebles! (l2004D:lDriver et all (120071): 
iMenard et all l|2010f) : iFukugita et all (|201lD . We have 
combined the points from Menard et al. (2010) for 
the dust contributions of halos and those from Menard 
& Fukugita, (2012) to a single set of data points, im- 
der the assumption that the amount of dust in halos 
doesn't evolve significantly with redshift. This is an 
assumption and could be tested in future data. At 
high re dshifts our estimate is co n sistent with the re - 
sults of iZukugita & Feeble^ (p00l : lD7wer et all (|2007f) : 
iMenard et al. (201(J). 

We note that the Menard et al. (2010, 2012) measure- 
ments assume a reddening law appropriate for the SMC. 
A Milky Way reddening law would have resulted in a 
a factor of 1.8 higher dust masses, and thus dust abun- 
dance, than the values shown in Fig. 12 (see discussion in 
Menard et al. 2012). Since the extinction measurements 
make use of a reddening law consistent with SMC while 
the direct emission measurements of the dust abundance 
that we show here assumed a dust mass absorption coef- 
ficient that is more consistent with the Milky Way, it is 
interesting to ask why the two measurements shown in 
Fig. 12 agree. The UV reddening law related to extinc- 
tion based dust abundance estimates comes from small 
grains that dominate the absorption and scattering sur- 
face area. SMC differs from other galaxies in that it 
does not show a prominent 2200rA feature, which is as- 
sumed to come from carbon bonds. On the other hand, 
the far-IR emission that we have detected is likely domi- 
nated by large grains, usually assumed to be a mixture of 
silicates and carbonaceous grains. The difference in the 
reddening law between SMC and Milky Way then should 
not complicate the abundance estimates since extinction 
and emission may be coming from different populations 
of dust grains (e.g., Li & Drainc (200_2,)). 

While Fig. 12 is showing that the dust abundances 
from extinction measurements are consistent with direct 
emission measure from far-IR background fluctuations, 
the above discussion may suggest that this comparison 
is incomplete. It could be that this agreement is merely 
a concidence of two different populations. Thus the to- 
tal abundance of the dust in the Universe is likely at 
most the total when summing up extinction and emission 
measurements. However, a direct summation of the two 
measurements is misleading and likely leads to an overes- 
timate. While small and large grains dominate extinction 
and emission, respectively, the two effects are not exclu- 
sive in terms of the different populations of dust grains. 
Some of the grains associated with extinction must also 
be responsible for emission. 

The far-IR background anisotropy measurements we 
have presented here have the advantage they capture 
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the full population of grains responsible for thermal dust 
emission in galaxies. The extinction measurements, how- 
ever, are biased to clean lines of sights where the lines 
of sights do not cross the galactic disks. We have cor- 
rected for the missing dust in disks by adding the density 
of dust in disks at z '--^ 0.3 to all measurements at high 
redshifts, but the disk dust density could easily evolve 
with redshift. The agreement we find here between the 
two different sets of measurements may, however, argue 
that there is no significant evolution in the dust density 
in galactic disks. In any case we suggest that one does 
not derive quick conclusions on the dust abundances or 
the agreements between extinction and emission mea- 
surements as shown in Fig. 12. There are built-in as- 
sumptions and biases between different sets of measure- 
ments and future studies must improve on the current 
analyses to understand the extent to which extinction 
and emission measurements can be used to obtain the 
total dust content of the universe. 

While the Herschel fluctuation measurements have the 
advantage we see total emission, they have the disadvan- 
tage that we cannot separate the dust in disks to diffuse 
dust in halos that should also be emitting at far-IR wave- 
lengths. In future, it may be possible to separate the two 
based on cross-correlation studies of far-IR fluctuations 
with galaxy catalogs and using stacking analysis, espe- 
cially for galaxy populations at low redshifts. These are 
some of the studies that we aim to explore with the H- 
ATLAS maps in upcoming papers. 



7. CONCLUSIONS 

We have analyzed the anisotropics of the cosmic far- 
infrared background in the GAMA-15 Herschel-ATLAS 
field using the SPIRE data in the 250, 350 and 500 //m 
bands. The power spectra are found to be consistent 
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with previous estimates, but with a higher amplitude 
of clustering at 250 /xm. We find this increase in the 
amplitude and the associated increase in the shot-noise 
to be coming from an increase in the surface density of 
low-redshift galaxies that peak at 250 /.im. The increase 
is also visible in terms of the bright source counts of the 
H- ATLAS GAMA fields (e.g., Ribgy et al. in prep). 

We have used a conditional luminosity function ap- 
proach to model the anisotropy power spectrum of the 
far-infrared background. In order to fit H- ATLAS power 
spectra at the three wavebands of SPIRE wc have 
adopted the spectral energy distribution of a modified 
black body and constrained the dust parameters Td and 
Pd using a joint fit to power spectra at 250, 350 and 
500 ^m. The results of our fit substantially confirm pre- 
vious results from the analysis of Herschel data and allow 
us to improve the constraints on the cosmic dust density 
that resides in the star forming galaxies responsible for 
the far-infrared background. We have found that the 
fraction of dust with respect to the total density of the 
Universe is ^dust ~ 10^^ to 8 x 10^®, consistent with 
estimations from observations of reddening of metal-line 
absorbers. 
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TABLE 2 

Angular power spectrum measurements at 250, 350 and 500 /im from GAMA-15 field of H- ATLAS. We tabulate the values 

AS PCi/2-n WITHOUT SHOT-NOISE SUBTRACTED. 
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